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Abstract 

We demonstrate that the partial wave decomposition of three-nucleon forces 

used up to now in momentum space has to be necessarily unstable at high 

partial waves. This does not affect the applications performed up to now, 

which were restricted to low partial waves. We present a new way to perform 

the partial wave decomposition which is free of that defect. This is exemplified 

for the most often used 27r-exchange Tuscon-Melbourne three-nucleon force. 

For the lower partial waves the results of the old method are reproduced. 
PACS numbers: 21.30.+y, 21.45. +v, 24.10.-i, 25.10,+s 
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I. INTRODUCTION 



Three-nucleon forces (3NF) act for more than two nucleons. The interesting question is 
their strengths and their signatures. 

A first observable where they clearly show up is the binding energy of the triton. Here 
it is known that the most recent realistic nucleon-nucleon (NN) forces cannot produce the 
triton binding energy and 3NF's are needed in order to get the experimental number JO 0. 

A next logical step is to look into the 3N continuum. Our results based on numerically 
precise solutions of the three-nucleon (3N) Faddeev equations and realistic NN forces agree 
overall very well with experimental data |§. In elastic Nd scattering there is only one dis- 
crepancy, which clearly sticks out, the low energy vector analyzing powers, which depend 
sensitively on the 3 Pj NN force components. There is either an ambiguity in their determi- 
nation from NN data or one really sees a 3NF effect. The inclusion of the 3NF's, however, 
which have been worked out up to now, do not diminish that discrepancy || |4| . Right now 
it remains a puzzle . At very low energies, near the Nd threshold, one has to expect to see 
3NF effects, which are linked to the accompanying shift of the triton binding energy. But 
we predict also scattering observables, which do not scale with the triton binding energy || 
||. This is an interesting energy region to be studied experimentally. At higher energies up 
to about 60 MeV the 3NF effects we find for cross sections are mostly small (of the order 
of a few percent) || f| and sometimes even shifting theory in the wrong direction. For 
certain spin observables there can be larger effects |§. We apparently stand at the very be- 
ginning in understanding 3NF's. At much higher energies, however, above 100 MeV nucleon 
lab energy, we predict big 3NF effects, though still based on a truncated partial wave basis 
for the 3NF 0. 

The present article is meant to provide the theoretical and numerical ground to add more 
partial waves to the 3NF. With the previous algorithm || for decomposing the 3NF into 
partial waves we ran into unexpected numerical instabilities for higher angular momentum 
states. 
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In section [TJ we show the origin of this defect using a model 3NF. We demonstrate that 
the failure in the partial wave decomposition (PWD) used up to now cannot be cured for 
higher angular momenta. 

In section [TTT] we introduce a new PWD. Thereby we decompose the 3NF into two quasi 
two-body parts and a recoup ling of Jacobi variables in between. This leads to a new PWD, 
which is free of the inherent difficulties of the old PWD. We perform the new PWD for 
the 27r-exchange Tucson-Melbourne (TM) 3NF [|K|] ||. Thereby we also comment on the 
numerics where it is relevant. 



We demonstrate in a numerical study in section [TV| that for lower partial waves the old 
and new PWD give the same results for Nd scattering. In a bound state calculation we 
extend the new PWD to higher partial waves and reproduce exactly the result of the Los 



Alamos - Iowa group |TT 



Finally we summarize in section [V} 

Appendix A contains the expressions for the permutation operators which are needed in 
our new approach. They are generalizations of the permutation operators which we usually 
use [[T2|] [[HJ in the Faddeev scheme. 

The permutation operators lead to interpolation of very many matrix elements. This 
includes a two-fold interpolation, as is shown in appendix A, for which a very efficient and 
accurate algorithm has to be used. It turned out, that our usual scheme, namely cubic 



basis splines ||TJ] , requires too much CPU time, which would not allow to calculate the 3NF 
within a reasonable time. Therefore we use cubic hermitian splines which need much less 
operations for one interpolation than the basis splines. In appendix B we briefly review the 
introduction of cubic hermitian splines and rewrite them to the most efficient form for our 
purposes. 
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II. THE PITFALL IN THE OLD PARTIAL WAVE DECOMPOSITION 



The 27r-exchange 3NF is shown in Fig. [l|. The three diagrams in Fig. [I] differ only by 



cyclical and anti-cyclical permutations of the three particles. Therefore it is sufficient to 
calculate the first diagram and find the full 3NF by the application of proper permutation 
operators. The first diagram is denoted by V} since particle 1 is singled out. 

We need the matrix elements of V4 in a basis \pqo), where p and q are Jacobi momenta 
(p is the relative momentum between particles 2 and 3 and q is the momentum of particle 
1 relative to the pair 2-3) and a denotes discrete quantum numbers. In detail 



where Z, s and j are orbital angular momentum, spin and total angular momentum in the 



denote the conserved total 3N angular momentum and its magnetic quantum number. In 
iso-spin space t refers to particles 2 and 3, which couples with 1/2 to total iso-spin T and 
its magnetic quantum number Mt- Thus the matrix elements we are looking for are of the 
form (pqa\ V} \p'q'a'). 

We found that for p = p' large, I — I' > 4, q — q' large and A = A' > 4 the matrix elements 
(pqa\ V} 1] \p'q'a') are some orders of magnitude larger than neighbouring non symmetric 
matrix elements. The increase of p = p', q — q f , I = I' or A = A' leads to an explosion 
in magnitude of the matrix elements. Moreover, we found that in those partial waves the 
symmetry of {pqce\ v} 1 ^ \p'q'a') under exchange of the primed and unprimed variables, which 
is manifest by the operator structure, is violated. This symmetry violation happens even 
for smaller values of Z or A. 

Such a behaviour hints clearly to a numerical instability of the code. This instability is 
not caused by an error in the code nor in the formal expressions of the PWD, but by the 
numerical realisation on a computer which always relies on a finite number of digits. We 
shall now illustrate the problem for a model using the old PWD. 




(1) 



two-body subsystem 2-3. A, \ and / are the corresponding quantities for particle 1, and JM 
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The instability is connected to the values of the momenta p, p', q and q' and the cor- 
responding angular momentum quantum numbers I, V, A and A' and is not related to spin 
or iso-spin. Thus we can use a spin and iso-spin independent 3NF model for our purpose. 
Moreover, the behaviour in q and A turns out to be analogical to p and I and therefore it 
suffices to treat only the p and / dependence in the 3NF model. (The old PWD for the full 
TM 3NF, which we used up to now, can be found in || and in more detail in fl5| .) 

Let us consider v} as depicted in Fig. |[ Neglecting spin and iso-spin, meson-nucleon 
form factors and replacing the shaded blop by a constant, we are simply left with the two 
pion propagators 

Vi 1] => 7^ ~~^~2~ 7 (2) 



Q 2 + ml Q> + ml 



where Q and Q' are given by 



Q = p-p> - -(q-q 1 ) 

Q' = p-p' + \{q-q') (3) 

We further simplify by freezing the ^-dependence and taking only one of the two prop- 
agators. Since the problems occur for p = p' we choose \p\ = \p'\ and allow only for an 
angular variation. Thus the model 3NF becomes 

ymodel 

(p(j) — p') + a) 2 + ml 

where a is a fixed vector. 
We define the two angles 

x 2 = p ■ p' ; X\ = (p — p') ■ a (5) 

and rewrite eq. (|) as 

ymodel ' 



2p 2 (l — x 2 ) + a 2 + 2apy 2(1 — x 2 )x\ + m 2 
1 A- 1 



B + Ax x f + xi 



(6) 



with 



A = 2a Py /2(l - x 2 ) 

B = 2p 2 (l-x 2 )+a 2 + ml 



(7) 



The expression in eq. (|6|) is expanded in Legendre polynomials resulting in 

oo / R \ 

Here and in the following we use the abbreviation 



(8) 



1 = 21 + 1 



(9) 



for orbital and spin angular momenta. 

The Legendre polynomials themselves can be rewritten in terms of coupled spherical 
harmonics with respect to the two directions involved in x\\ 

Air 



(10) 



with 



(ii) 



Further one has 



Yr(a+b)= y: 



i Xl b X2 



x 1+ x 2 =i \a + b\ l \ A x ! A 2 ! 
We insert now eq. ([U]) and eq. ( |T2"D into eq. @ and get 

'B\ Air 



Air V i 



(12) 



3C 



model 4—1 



A' 1 E h Qh ( -j 



p Xi (-p'y 



x E 



4tt /i! 
\ Ai! A 2 ! 



{yt^P%Y h {a)} 



00 



(13) 



Hereby we used the notation 



{a 2l , 6/ 2 }' m = C (hv km - u, Im) a v h b 



m—v 
h 



(14) 
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to describe the coupling of two quantities a and h with angular momenta l\ and l 2 to Im. 
Already here we can see the problem with this type of PWD. The expression in eq. ( |13|) 
seems to have a singularity for p = p', though the model force is not singular, of course. 
The singularity in eq. (^) is only apparent and is removed by the sum over Ai + A2, if 
that sum can be evaluated accurately. The next step in the PWD is to expand the x 2 - 
dependence in terms of Legendre polynomials analogically to what we have done for the 
xi-dependence. The x 2 dependence occurs in the denominator \p — p'\ [l and in the Legendre 
polynomial of the second kind, Qi(B/A). Performing that expansion we separate the singular 
denominator from the sum over Ai + A2, which has to cancel the singularity. This separation 
and the necessarily limited accuracy in the numerical evaluation of the sum over complicated 
geometrical expressions is the reason for the numerical instability as we shall demonstrate 
now. 

Let us now perform the x 2 expansion: 



1 00 I — Att 1 I 

vr del = i- e e R\ fHn {y'Ufofaw)} 

i 1= o \i+\ 2 =h N A i ! A 2 ! 

00 , — 

x £ 2^i 2 {-r VZJPJ) H hl2 {p,a) (15) 
The expansion coefficients Hi^ip^a) are given by 



2p' 2 (l-X2)+a 2 +ml 



H h i 2 (p,a) = I ' dx 2 P h (x 2 ) — - v ; 2ap ^ 2{1 / (16) 



,+1 Qh 

dxo PiAxo) — 

(^271^X2 

The rest of the PWD is straightforward: one has to recouple the spherical harmonics in 
eq. ( |15D in such a way, that Yx^p) gets coupled with Y i2 (p) and Y\ 2 (p') with Yi 2 (p'). The 
intermediate form can be simplified using 



= c(ho / 2 o,/o) Yr(d) (17) 

Finally projecting onto definite orbital angular momenta one ends up with 



(plm\ Vr dd \v'l'm!) = tt E yh E V4 H hl2 (p, a) Yf™' {t 

a P h h 



4"7T h\ 



Z 2 Ai A 2 



h Ai A 



2 



X / a — - y t 2 /\]y\ 2 \ 

x 1+ x 2=h \ Ai! A 2 ! [ j 2 r / 

x C(Z 2 AiO,ZO) C(/ 2 A 2 0,/'0) 

x C(hm' h - mf, 00) C(lm I'm' - m, hmf) (18) 

We see from eq. ([IB] ) that the p-dependence is now shifted into the expansion coefficients 
f/^ 2 (p, a). Since the code for the 3NF failed for large p we regard the limes of (p, a) 
for p — > oo. Inspection of eq. (|T6| ) suggests to introduce as a new integration variable the 
argument of : 

z( X2 ) = 2p2(1 ~" 2)+a2 + m - (19) 
2<W2(1 - x 2 ) 



Because z{x%) has a minimum at x 2 = 1 — a the integral over z splits into two parts. It 
results 

. /„2_i_„,2 / 

2a 2 z a 2a 2 z 2 — (a 2 + m 2 * 



H hh (p,a) = / 4p2+a a 2+m| tfo 



4ap 



P 2 ^a 2 ^ 2 - (a 2 + m 2 ), 



a 2 + m 2 a 2 z 2 az 



a 2 + m 2 a 2 z 2 az 



ii+i 

2 



x ^) (2 1^-^^ + -r + -f V^ 2 - (« 2 + <) J J 

/•°° , / 2a 2 z a 2a 2 z 2 - (a 2 + m 2 ) \ 

a 2 + m 2 a 2 z 2 



„ / a~ -+■ m: a~z~ rr - ; 77 77 



t]+l 
2 



(21) 



Because Qi x {z) decreases like z~ lx ~ x for large z, z is limited to small values and both argu- 
ments of P[ 2 go to 1 for p — > oo. As a consequence f/^ 2 (p, a) becomes independent of / 2 . 
We end up with 

lim H hh (p,a) (xp 11 - 1 (22) 

p— >oo 



Because the model force decreases for large momenta the increase of H for p — > oo has 
to be canceled by the sums occurring in eq. (|i"8"|). Indeed this is the case, as we will show 
now. 

If H is independent of Z 2 we can perform the sum over Z 2 analytically: 

. f h A x A 2 
V/ 2 <^ } C(Z 2 AiO, ZO) C(l 2 A 2 0,/'0) 

; 2 U V I 



(-) 



A 2 



v - C(AiO A 2 0,/i0) C(hO Z'O, ZO) 
\ h 



(23) 



Using this result we can also perform the sum over Ai + A 2 analytically: 



E ( 

Ai+A 2 =ii 

E < 

Ai+A 2 =ii 



\A 2 



\ A 2 



4tt (2/i)! 
\| (2A X )! (2A 2 )! 



C(AiO A 2 0,/i0) 



4vr (2ZQ! 
\J (2A0! (2A 2 )! V 1 \ 



(2Zi-2Ai)! (2Zi-2A 2 )! 



(2k + 1)! 



(Zi - Ax)! (/i-A 2 )! 



47T 

Ai 



Ai / 1 (/ 1 -i)...(/ 1 _A 1 ) 



Ai! 

= (-)' 1 v 7 ^ (1 - l)' 1 



(24) 



h 

\i \Ai 

From this result we can see two things: firstly the term (1 — l)' 1 cancels indeed p' 1 ^ 1 . (A 
somewhat harder inspection shows that only p~ 2 survives, as it should.) But secondly, to 
achieve this cancellation one has to add up an alternating sum. This causes loss of digits if 
the consecutive numbers are of similar magnitude. This is the case here. For large l\ it will 
be impossible to get a reasonable result for (1 — l)' 1 . 

Looking into the corresponding sums in the code for the full TM 3NF we found, that 
for diagonal matrix elements with p > 6 fm _1 , I = 4, q > 12 fm _1 and A > 4 this cancel- 
lation occurs with consecutive numbers being of alternating sign and equal within up to 8 
digits. Thereby we lost all numerical accuracy. As a consequence the increase from 
the corresponding H could not be compensated and badly wrong results appeared. 

It is also clear that an increase of the numerical accuracy cannot solve that problem. 
Higher numerical accuracy would only shift the problem to higher I and A. 
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III. THE NEW PARTIAL WAVE DECOMPOSITION 



In this section we will present a new PWD, which avoids the introduction of artificial 
singularities. The new PWD leads to a numerically stable code. We perform this PWD for 
the TM 27r-exchange 3NF. In the notation from Fig. this force is given by P 



X 



(2tt)6 Ami Q 2 + ml Q' + m 2 

t 2 ■ r 3 (a + b Q ■ Q' + c (Q 2 + Q' 2 )) + d n • f 2 x f 3 a x ■ Q x Q> 



(25) 



The constants a, b, c and d are determined by low energy theorems [|T(J]. The form factors 
are parametrized as 

F(Q 2 ) = ^< (26) 
A 2 + Q 2 

with the cut-off parameter A. 



The Faddeev scheme [|l3j \12\ JT6| JT^] we are using requires matrix elements 



(pqa \ Vj 1 ^ (1 + P) \p'q'ce'), where P is the sum of a cyclical and an anti-cyclical permutation 
of the three nucleons. Henceforth we add a subscript 1, 2 or 3 to the basis states \pqa) in 
order to distinguish the three different choices of Jacobi momenta, which single out either 
particles 1, 2 or 3 to be the "spectator nucleon". In that notation we can write 

.(pqalV^ (1 + P) Ip'q'a'), 
= Moc\ V} 1] \p"q"a") 3 3 (p"q"a"\ (1 + P) Ip'q'a'), (27) 

We introduced intermediate states of the type 3, since they single out the pair (12), which 
participates in one of the two pion exchanges in the operator V} . This is the natural 
basis to evaluate that specific pion exchange. The other pion exchange is in the pair (13), 
for which the basis \pqa) 2 is the natural one. This leads us to the basis idea of the new 
PWD, to split each term of V4 into two parts, which are related to the separate one pion 
exchanges and which are then necessarily linked by a change of coordinates. We call the 
corresponding two parts f| I3 and u| I2, where the subscripts stand again for the spectator 
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nucleon. Thus each of the four parts of V4 in eq. (|25|) , called in the following, can be 
cast into the form 

^pqa\ v ( l ] \p'q'a'} 3 = ^ ^ ^ x {pqa \ p"q"a") 2 
x 2 {p"q"a"\ vj£\ \p'"q"W") 2 2 {p"'q"'a"' \ p""q""a"") 3 
x 3 (p""q""a""\v£ ) \ 3 \p'q'a') 3 (28) 

The advantage of this PWD is that the matrix elements 2 (p"q"ct"\ U4 I2 \p"'q"'a"') 2 and 
(p""q""a""\ I3 |pV a ')3 are easily evaluated, and neither these matrix elements nor the 
two changes of basis states introduce numerical instabilities. These changes of basis states 
are performed by the application of certain permutation operators, which are given in ap- 
pendix A. Their actual use requires 2-fold interpolations which are performed with the help 
of cubic hermitian splines (see appendix B). 

The basis states \pqa) 1 in eq. (^7]) are antisymmetric in the two-body subsystem (23). 
This is required in the context of the Faddeev equations |L3| [0 |16|] The states \pqa) 2 



occurring in eq. (^) do not have that property, since the permutation operators remove 
that symmetry and therefore one has to sum without restriction over all quantum numbers 
to a given total angular momentum J and parity. On the other hand the states of type 
3 have again to be chosen antisymmetrical in the corresponding two-body subsystem (12). 
This follows from eq. (|27|). The permutation operator (1 + P) applied onto \p'q'ot') l in 
eq. ( |2"TD generates a totally antisymmetrical state and consequently the intermediate states 
3 (p"q"a"\ have to be antisymmetrical in the subsystem (12). The expression (^) can further 
be simplified. One has 

Pis P2s\pqa) 1 = \pqa) 3 (29) 

and 

P 23 P 13 (1 + P) = P 23 P 13 (1 + P 12 P 23 + P 13 P 23 ) = 1 + P (30) 

Consequently 
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.(pqalV^ (1 + P) \pWa') x 

f, 1 (pqa\Vj 1) \p"q"a") 3 x {p"q"a"\ (1 + P) \p'qW) 1 



(31) 



The basis states defined in eq. ([I]) are a direct product of space-spin states and iso-spin 
states: 



\pqa) = \pqaj) \a T ) 



p{ls)j q(X-)I (jI)JM 



(t-)TM T 



(32) 



The iso-spin dependence of v± factorizes off and one is therefore naturally lead not to 
include the iso-spin dependence into the splitting of and to evaluate the iso-spin matrix 
elements 2 (a T \ . . . \a' T ) 3 directly. One easily finds 



2X it-)TM T 



t 2 • r 3 



{t'\)T'M' T \ = 5 TT , 5 Mt m t , (-6) (-)* \Jti' { 

* I 3 



1 1 ,i \ 

2 2 1 

1 1 ± 

2 2 

t \ T 



(33) 



and 



2 ,(t-)TM T 



ti-t 2 x r 3 



it'\)T'M' T/ 3 



5 XT' o 



M T M T i 



24i (- 



>2T 



it' E ! 



A ± 1 



1 1 

2 2 



t 



' T 


i 

2 


t ' 


< \ 


1 


A » 


t' 

V 


1 
2 


i 

2 J 



(34) 



These matrix elements will be denoted by 2 ( a r| I a WtIz in following. Restricting 
the splitting of then to the space-spin parts and calling the two parts u>| | 2 an d w^' | 3 
eq. (p8|) achieves the form 

x ( 2 (P // 5 //Qt jl w 4 b \p\q\OL\j) 2 2 {piqia 1J I p 2 q2a2.j} 3 
x 3 (P292Q!2j|wl 1) |3|pVa'j)3 ) 

x 2 «| J« \a' T ) 3 (35) 

In eq. ( |35| ) we introduced intermediate sums ^ and where the related quatities carry 
the indices 1 and 2. Since the quasi two-body forces W4 |2,3 are diagonal in the spectator 
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quantum numbers, the intermediate sums over a±j and a 2 j are quite restricted and the 
main introduction of auxiliary states in eq. (|35|) are the ones with a", which moreover are 
unconstrained by any symmetry. The a' states are anti-symmetrical as mentioned above 



but otherwise unconstrained due to eq. fl3~T|). We should also note that the operators w^\2. 



3 



can change the parity and the total angular momenta by ±1, as will be shown below. 

The evaluation of the permutation operators in eqs. (|3l|) and (|35|) consumes essentially 
all the CPU time; the time for the calculation of the w^\ 2 ,3 matrix elements is negligible. 

Let us now regard the matrix elements for w^\ 2 ,3 related to the four terms proportional 



to the constants a, 6, c and d in eq. (^). We omit the common factor ^— ^ an d the 
constants a, b, c and d. 

The a-term naturally splits into two parts, one related to a pion exchange between 
particles 1 and 3, 

W W\« F^Q*) ? 3 '® (36) 
Q 2 + ml 



and the other to a pion exchange between particles 1 and 2 

Q' 2 + ml 



^ ] \t = F\Q' 2 )4f^- (37) 



Consequently one has to evaluate 



M J / = 2 (pqaj\ F\Q 2 ) |p 1 g 1 a 1J ) 2 

Q 2 + m 2 

M J / = 3 (p 2 q 2 a 2J \ F\Q' 2 ) _* 2 ' Q ' \p'qWj) 3 (38) 

Q +ml 



The states atj from eq. ([32] ) couple the total angular momenta of the two-body subsystem 
and the third particle: 

r ^ \^ JM 

M = {W 3 )Wi)Y M ^{\{ls)3) (A-)/)] (39) 



Further the pion momentum Q from Fig. |2| is expressed in terms of the Jacobi momenta JT3 
in the (13) subsystem as 
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Q = p-p l 



(40) 



Apparently M^ 2 is diagonal in the quantum numbers of the spectator nucleon, and we find 
r j,2 _ S(q-qi] 



S X x 1 6 Ih ]T c U m IM ~ m ' JM ) c (h™i hMt - mi, JM) 



x 2 (pj m \ 



&3 ■ (P-Pl) 



(41) 



(p-PiY + ml 

The scalar 3% ■ (p — pi) formed out of two tensor operators of rank 1 can be decomposed 
using standard techniques |13 



as 



<?3 • (P-Pl) = -V3 {(7(3), (P-Pl)!} 00 

= -4vr £ P a (-Pi) b M3),?/; fe (p,Pi)} 



(42) 



a+b=l 



As in eq. ( 14]) the curly brackets denote angular momentum coupling and the product of 
two spherical harmonics coupled to angular momentum 1. Note that in contrast to eq. ( |T3"| ) 
there arises no inverse power of \p — pi\, which would lead to an artificial singularity with 
the fatal numerical consequences described in section [TJ. 

Next we expand the pion propagator and form factors into Legendre polynomials de- 
pending on x = p ■ p\ . We rewrite 

\ 2 



i 



A 2 - mi 



(p-Pi) 2 + m l \A- 2 + (p-Pi) 



(43) 



as 



with 



(A 2 ~ ml) 

{^PPiY 



i 



(Ba - B, 



a i * ) V B<m-K •' B\ X 



B\ — B mn (B A -x) 2 



(44) 



Br, 



p 2 +p\ + A 2 

Zppi 
p 2 +pl + ml 

2ppi 



(45) 



It results 
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OO 7 

1=0 



with 



i A 2 — m 2 

HfaPi) = — (Ql(B m J - Qj(B A )) + — ^ QL(B A ) 

These steps lead to the intermediate result for the matrix element in eq. (jO) 



M 3 a = 2 {pjm\ 
(4tt) s 



^3 • (P-Pl) 

(p — pi) 2 + m 2 



F 2 ({p-Pi) 2 ) \pijimi)r 



2 (p 3 m\ y: p a (-Pif w(3),yi b (p,pi)r 



a+b=l 



1=0 



i Hj(p,pi) y™(p>pi)\piji m i)-< 



(46) 



(47) 



(48) 
(49) 



or 



= -2nVS P a Pi E lH MPi) £(- 



V'2 



a+b=l / = g n«2 

a b 1 

a < !> C(aO Z0,ii0) C(60 Z0,i 2 0) 



%2 H I 



(MMsiylMp!)} 00 \j imi ). 



(50) 



This can be expressed in terms of reduced matrix elements: 

^5 



(jm\{a(3),y l i .(p,p 1 )} 00 \ 3imi ) 2 =S sh 5, 



mmi 



X 



ll Si j 

( _ 1} 1 +Wsa J 

V3 ( s I 1 

y l u Ap,Pi)\\k) 2 2 ( s |k(3)i| Sl } 2 



(51) 



with 



and 



,(s||a(3) ||si) 2 = V6 Js\ (- 



\s+l 



! 1 Si 

2 2 S1 



ytiMPl) I l Z l> 2 = 5 n« <w 



1 s 



NT 



(52) 



(53) 



In view of later use the expression ( p3| ) is formulated for general L. Note we define the 
reduced matrix element as 
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C(Lm - mi hmx, Im) {1\\O l \\k) = (lm\ O lm \hrrn) (54) 
The result for Mt 2 is 



Mf = 6tt V2 6 jh 5 mmi J§f x 



Is I J Is Z 1 



m ; f a 6 1 1 

x 2 P B Pi £i H j(p, Pl ) I _\ C(aO 10,10) C(bO 10, 1$) (55) 

a+f>=l 7=0 [h I I ) 

Finally we put this expression into eq. ([11]) for M^ 2 . This allows the summation over 
the magnetic quantum numbers m and mi. We also can perform the sums over a + b = 1 
and 7. We end up with 

Ma' 2 = ~ ' fijji ^AAi ^JJj $JJ 1 ^MM X <fy-ii|,l 

I I Si) (h *1 J 
IS i J 1 S / 1 



x2n V6 (-)J+ 1+S+S1 Jssi 



2 



x ^max(Mi) (p ^(p,pi) -pi fl,(p jPl )) (-r^ (56) 

This result for M^' 2 is a very simple expression, in contrast to the a-term using the old PWD 
(see P| and [IB]). The complicated structure of the 3NF in our new PWD is hidden in the 



three permutation operators of eqs. (]5ID and ([55]). The latter ones are, however, well under 
control algebraically (see appendix A) and numerically. The M^' 3 can be gained out of M^' 2 
by carefully replacing quantities with the result 



M J /{ V a, Pl ax) = (-)!+'+* M J /( Pl a uP a) (57) 

The splitting of the operator introduces parity changes for the individual terms M^ 2 
and M/' 3 , as is obvious from eq. (|56|). Their product, however, conserves parity, as it should. 



The c-term is very similar to the a-term, only Q 2 and Q' 2 factors occur additionally: 

<?3 • Q 

Q2 _|_ m 2 



M J / = 2 (pqaj\ Q 2 F 2 (Q 2 ) °'*' Q - \ Piqi a XJ ) 2 



M J / = 3 ( P2 q 2 a 2J \ Q> F 2 (Q> ) _" 2 2 '- V bYa'^ (58) 



2 >'-\ 

g /2 + m 2 

This leads to a replacement of H by 
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~ 771 ^ — TH^ 

Hj(p, Pl ) = (Q T (B m J - Q T (B A )) - _— * A 2 QL{B A ) (59) 

PPi 2{PPv 

in eq. (|56| ) and provides M^ 2 and M/' 3 . The symmetry relation eq. (|57|) is also valid for 

M C J ' 2 and M C J ' 3 . 

— * — * 

The b-term in eq. (|25|) has the scalar product Q ■ Q' on top of the structure of the a-term. 
That scalar product Q ■ Q' can be rewritten as: 

-> -» 4-7T 

q.qi = -^3 {Q u g;} 00 = -V3 -f qq' y JJ(g, go 

4-7T 

= gg' E(-) M Y i(Q) y?W) ( 6 °) 

This leads to the amplitudes 

M b J ' 2 = M 2 (pqaj\ Q Yf(Q)F 2 (Q 2 ) J^L \ Piqi a u ) 2 

M b J ' 3 = M 3 (p 2 g 2 « 2J | Q' Y^(Q')F 2 (Q' 2 ) % ' ® Ip'q'a'j), (61) 
V 3 Q> + m 2 



Now similarly to eq.([|l]) one has 



M J,2 = 5(q _gi) ^ ^ ^ JM _ ^ JM) c , ( ^ TOi /iMi _ m ^ JiMi) 

9 mmi 

x V"? 2 (kH 5 3 ^((p-Pi) 2 ) |p-Pi| Yfip-pi) \pijim x ) 2 (62) 

Its evaluation requires the following recoupling: 



Yfip-P!) 5(3) ■ = -v 7 ^ |p-pi| {^(p-pUM 3 ), ^-pi)}°} 



171 Ip-pii E(-r +1 Hsj^scr-A^-A)} 1 * 



3 

Hp-pII £(-) il+1 c(io io,iiO) {^^(p^)} 1 '" 

= \p-pi\ H3),w^o} 1,M -yf Hs),^^^} 1 '" 

As an intermediate step using eq. ( |T2"D we get 

' 2 - ■ It* l p j m \ ^3 • (P-Pl) 

3 2 (p-pi) 2 + m 2 



(63) 



M b 2 = \l— 2 {Pj m \ — — , f2 (p-Pi) 2 \P~Pi\ Y i(p-Pi) 
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— 2 {pjm (p - Pl ) — — — — — - 

3 [p-piY + m^ 

3 {p-pi) z + m^ 



W 3 )* S/SoCp-pi)} 1 ''' IP^mi) 



x 



E 



p a (-Pi) f 



W 3 )> 1/ ab(^~ Pl)} 1 ^ IPlJl^l) 



(64) 



The expansion into Legendre polynomials for the first term is the same as for the c-term 
and for the second term as for the a-term. We get 

r 2 



Mi = 



57T 



8n< 



E(-) T ^ Hj(p, Pl ) 2 (jm\ ygfah) {a(3), t/^^px)} 1 " 1 \j imi ), 
i 

AV15 J2(-Y ^Hj(p, Pl ) 



x 



E 



p a (—pi) 1 



a +b=2 

To proceed further we have to recouple 



(jm\ y™(p,pi) {<t(3), 2/^(p,Pi)} 1,A4 b'i"ii 



1>M 



i r-: . f 6 a 2 

= ^ VaME(-) l4+a+i . _ 

x C(aO 70, i 3 0) C(bO Z0,i 4 0) {(7(3), f/^Pi)} 



(65) 



(66) 



and 



2/™(P.Pi) H3).y«(P.Pi)r = i H3),2/J(p,Pi)} 



1,M 



47T 

which leads directly to the following two matrix elements: 



i(jm\ {(j(3), yl Ai4 (j>,Pi)} 



5 i3 i 5 iih C(lfi iimi, jm) 3^ (-) s m Vj'issi 



lie 

2 2 61 

i . i 



and 



f 2 1 1 

{ 1 s j J 



(67) 



(68) 



,(jm| {(j(3), y^p.Pi)} 1 ^ Ijimi) 



l l 



1 s 



j jl 1 

Si s I 



(69) 



and then to 

,2 



M 3 h ' = (7(1// jxmi, jm) 



x 



2tt i 



X 



40tt V6 Vj'isS! 



^lHj(p, Pl ) £ 



1 I s 

2 2 

1 a i 



lie 
2 2 *1 

1 S i 

r 2 i n 

il Si ji 



Si s Z 



J s 3 . 
b a 2 

C(aO Z0,Z0) <7(60 Z0,/i0) 



(70) 



a+b=2 yj(2a)\ (26)! [/ / x I 
The summation over the magnetic quantum numbers in eq. ( |62T) using the m, mi dependence 



from eq. (|f0|) can be done analytically: 

C(jm IM - m, JM) C{j x mi I x Mi - m x , J\M X ) (7(1// j x m u jm) 

1 Ji J 



_\l+jl+I+J 



(") 



/jJi C(l// JiMi, JM) 



I J Ji 



(71) 



J,2 



and we get for M'i 
r j,2 S{q-qi) 



Mi 



Q 



5 XXl 5 Ih CilfiJ.M^JM) 



x (-) 



i+j 



jjishJi 





1 

2 








{ 


S 




j 





X 



2tt 



^ ^ (-y +i fli(p jWL ) ' ; 11 1 

6 1 Si s I 



- 40tt v 7 ^ (-; 



s+ji 



r 2 i n 

Zl Si ji 

L * s i J 



6 a 2 



I a ^ 2 (26)! [/ / x I 

rJ,3 



C(a0 10,10) C(b0 10, kO) 



(72) 



The evaluation of M b ' proceeds along similar lines or can be abbreviated by judicious 
replacements and using symmetries of 6j- and 9j-symbols. The result is 
r J,3 S(q'-q 2 ) 



Mi 



S yX2 S rh C(l-fiJ'M',J 2 M 2 ) 
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x (") 



r+j 2 +s'-s 2 



1 1 e 

2 2 S2 



1 S 



/ 1 



i h f 
v J' j 2 



3 \s 2 s' I' 



40vr V6 (-Y'+i* 



r 2 i n 

^2 S 2 J2 

Z' s' f 



b a 2 



a+ b=2 s/{2a}\ (26)! [/' Z 2 I 
Now, looking back to eq. fl60D , we can perform the sum over \i: 



C(aO ZO, Z'O) C(60 ZO, Z 2 0) 



(73) 



C(l/i Jxmi, JM) C(l -fi J'M', J x Mx) = <W (-) 



Ji-j 



\1 



(74) 



Thereby we used J\ = J 2 and Mi = M 2 , which results from the fact that the permutation 
operator standing between M b J ' 2 and M 6 J ' 3 (see eq. (p5|)) conserves J and M. Thus despite 
the fact that M b J ' 2 and M fe J ' 3 do not conserve J, as is obvious from the Clebsch-Gordan 
coefficients in eq. fl74|), the total expression f| |& does it, as expected. The parity, however, 

J 2 .7 3 

is conserved by the individual quantities M b ' and M b ' . 

Having performed the /i-summation reduced M-quantities occur, which are defined as 



M/' 2 = C(lfi JiMi, JM) (-) Jl ~ J M 



.1.2 
b 



M fe J ' 3 = C(l-/i J'M',J 2 M 2 ) 
and which obey the following symmetry relation 



u^ 3 



(75) 



j 2 

M 6 ' (pa,piai) 



Sl — s 



F J,3 



^b' (piQ!i,pQ;) 



(76) 



Again we see that the expressions using this new PWD, now for the 6-term, are much 
simpler than the expressions of the old PWD, see and ||15|| . 

The last term, the cZ-term, contains the operator a\ • Q x Q' (see eq. (p5|) ): 
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= -iV2^QQX(-)" {^(l)^i(Q)} 1 '^ Yf(Q') 
= zv / 2^gQ'E(-r Y^')}^ yfiQ) (77) 

According to the two possibilities to couple cr(l) we have two ways to write down the d-term: 

l (vq*\v { l ) \ d \p'q'a'), = P^ 2 M J / P 2 ^ M J / 

= Pi~2 E(") M M b' 2 P2~s Mf (78) 

with 



Mf = -iV2J^- 2 (pqaj\Q UllY^Q)} 1 " F 2 (Q 2 ) | Ml(llJ ) ; 



^v^/- 3 (p 2 g 2 a 2J |Q' {(7(1), nC*)} F 2 (Q' ) J?'* [pV^) 3 (79) 
1 1 1 J Q +ml 



The permutation operators P±^2 and P 2 ^3 in eq. ( |fq) stand for the recoupling matrix 
elements (see eq. (|35f)) and we dropped the iso-spin factor I4 . We see that M^' 2 and M^' 3 
are purely imaginary. Together with the imaginary iso-spin factor of eq. (PH) , the matrix 
element is of course real. 

For Mj' z we have to recouple the operator 

a 3 -Q {<T(l),Y 1 {Q)} 1 * = # 3 -(p-i> 1 ) { ( r(l),n(p^p 1 ) V "' 



^3 IP-Pil {{tT(3),y 1 (p-p 1 )}°,{(7(l),F 1 (p-p 1 )} 

z 2 ii 1 



1 1 1 



= |^-pi| E(-) i2 Vii C(10 10,i 2 0) 
x {Ml),a(3)}\y ia (p^p 1 )} 1/ ' 

= -L|^-^l {mi),^)} 1 ,^^)} 
+ V2 Ev^r 11 ) {Mil.^.y^M)} 1 " (so) 

H = l,2 [ 1 1 1 J J 

In the last step we inserted explicitly the sum over the two i 2 -values. Finally we decompose 
the spherical harmonics and get 
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a 3 -Q {(7(1), Y 1 (Q)} 1 " 

= ^ \p-pi\ {hi),^)} 1 , yl&Pi)} 1 " 

, 1 pr f2 h 1 

+ 8v / 15^ ,^ ^ V Vii ^ 

b-Pil ut^, 2 I 1 1 1 



E {wi)^(3)r, 2/Lcp.pi)} 1 " ( 81 ) 



X 

a+b=2 \la\b\ 



7 9 

which inserted into M d ' yields: 

Mf = -iV2 f 2 {p qaj \ ip- p.i 2 {wi),^)} 1 , y ^(p.pi)} 1 " ipi«i« u > ; 



il=l,2 [1 1 1 J a+6=2 y a!6! 

X 2(Pgaj| |^-^+ 1 m^ {Wll.^f^K^OriM^ (82) 

This expression contains no singularity. 

It remains to expand the ^-dependence in terms of Legendre polynomials, which leads 

to 

Mf = -iV2^ J2\/t (-) T H T (p, Pl ) 

6 I 



X 2 



2/J(p,pi) {{(t(1),(7(3)} 1 , J (]3,pi)} lM Ipigiau) 



i327r 2 VlO 



E A ( 2 " 1 ) E ^p=- E ^ #t(p.pi) 

U=l,2 [1 1 1 J a+fe=2 A/a!6! I 



x 2 (m«j| Vl&Pi) {W( l ),v(3)} n , ylbiPiPi)} M bigi«u) 2 (83) 
and to combine the spherical harmonics 

t/gWi) {Hi),^)}-,?/^,^)} 1 " 



{Wl)."(3)} il ) |/i(p 1 Pi)} 1 2/s(p ) ft) 



1/i 



= V (_) I + J 4+« C(aO 70,i 3 0) C(bO I0,i 4 0) 

x 1° 6 2 | (Hl),^)}-,?/^^^)} 1 " (84) 
[U H I J 
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and 



t/S[(p,pi) {WW, vmViyl&pi)} 1 " =^ {wi)^(3)} 1 ,2/j(p,pi)} 1,i 

Thus we have to evaluate two matrix elements. The first is 
2 (pqa\ {{a(l),a(3)} H .y^p.px))^ bigi«i) 2 



mmi 



x Cilfijmnjm) (-) h+1 V^S < 



r 2 ix 11 

Jl Si jl 



I S J 

x2(i\\ylMpi)\\h) 2 2 (siiwi),a(3)nis 1 ) 2 



Only the reduced spin matrix element is new: 

2 <s||Ml),a(3)n|si> 2 = 6^ < 
After summation over m and mi we get 



f 1 1 h\ 

lie 

2 2 

V 2 2 S 



2 {pqa\ {{(j(l),c7(3)} n , y l 3i4 (p,pi)} 1 ' 1 \piqiai) 
^ q 9 ^ 5 AAl S Ih 8 i3 i 8 iih 6V 



x C(\n JiMi, JM) 



1 ji 3 
I J Ji 





-i+i- 


-J+h ^ 


jjlSSiJih 


' 2 


ii 


1 ' 




' 1 1 h ' 


< h 


Sl 


jl 


> < 






s 


J . 




i 1 S 
I 2 2 6 J 



The second matrix element is 



2 (pH {{o r (i),^(3)} 1 , yl&Pi)} M bigi«i) 



S+I+J yjjissiJi 



x C(l// JiM u JM) 



1 Ji J 
I J Ji 







' l 


1 


1 ' 


n 




i 

2 


1 
2 






\- 


Sl » 






1 

^ 2 


1 

2 


s , 
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We end up with 



M, 



J,2 _ ^(g - gl) r r / N/+J+1 

" -3 "AAi d Ih {-) 



JJishJi C(l/i J X M X ,JM) 



1 3x 3 

I J Jl 



X 



5 Ml i 47r \/6 (- 



\l+s 



I S j 
1 jl Si 



r 1 1 n 

1 1 



2 2 Sl 
1 1 



> Hi(p,pi) 



+ i 240tt V6 (-) J1 J2(~) h *i 



2 2 

2 ?! 1 

1 1 1 

a 6 2 



' 2 


ii 


1 ' 




' 1 


1 




< h 




Ji 


> < 


1 

2 


1 
2 


Sl > 




s 






1 

*. 2 


1 
2 





x V . = yi yiH 7 (p :Pl ) { } C(aO 10,10) C(W 10, hO) 

a+b=2 y/(2a)\ (26)! j [h I V 

Similar steps or using judicious replacements lead to 
j,3 S(q'-q 2 ) 



M. 



Q 
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v/'+J 2 +2'+/2+l .hi 



X 



j'32s's 2 J' C(l -// J'M', J 2 M 2 

r 1 1 n 

( V s' ?' 1 



1 32 s 2 



1 1 c 

2 2 52 



1 1 

2 2 



1 J2 f 

V J' J 2 
Hi.(p'p 2 ) 



+ i 240tt V6 (-) h Y,( 



2 ii 1 
1 1 1 



' 2 


«1 


1 ' 




' 1 


1 




< Z 2 


«2 


32 


> < 


1 

2 


1 
2 




V 

V 


s' 


3' . 




1 

*. 2 


1 
2 


s' , 



E 



p/o pb 



a b 2 

Z 2 /' I 



C(aO 10, 1'O) C(bO 10, l 2 0) 



a+b =2^(2a)\(2b)\ j 

The last step is to perform the sum over /x, eq. (]7q) . Because in eq. 
together with M&, we have to define consistently with M& (eq. (|75|)): 



(90) 



(91) 



Md is used 



Mi' 2 = C(l/i JxMi, JM) (-) Jl_J 



M, 



and we have the relation 



d 

J,3 



C(l -/i J'M',J 2 M 2 ) 



\ J 2 



(92) 
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Mf (pa, Pl a 1 ) = (-) l+h Mf(p 1 a 1 ,pa) (93) 



Then the summation over /i can be done analytically using again eq. (74). 

This concludes the presentation of the new PWD for all four terms of the TM 3NF. In 
this new approach one avoids the pitfall of the old PWD and one can evaluate now the 3NF 
for higher partial waves. 

IV. NUMERICAL EXAMPLES 

This section has two aims. The first is to show that the old and new PWD give the same 
results for partial waves where the old PWD is still numerically valid. The second is to show 
that the extension of the new PWD to higher partial waves gives correct results. Lower 
and higher partial waves refer to the outer states of type 1 in the matrix elements of the 
3NF. The number of intermediate states of type 2 or 3 in eqs.([3Tl) and ( |35"1) are principally 
unlimited and the necessary number of states for an accurate representation of the 3NF 
matrix elements have to be determined numerically. We find, that in eqs. fl3T|) and (^) the 
a" sums require two-body angular momenta up to j max = 5. As already mentioned the 
additional sums in eq. (|35| ) over ay and a^j are essentially trivial, since the corresponding 



quantum numbers are only marginally changed by the action of w\ 1 2,3 from the ones of the 
states of type 2 and 3. 

We calculate the 3NF at 16 p-points below 10 fm" 1 and 16 g-points below 20 fm -1 . 
According to our experience this is enough to describe the 3NF [T2|] |T3|. The inner basis 



states of type 2 and 3 require 20 p- and g-points in order to achieve a fully converged result. 
The p-range is extended to 16 fm -1 in order to avoid extrapolation. 

As mentioned in section |1| the old PWD fails for I > 4 and A > 4. Therefore the largest 
two-body angular momentum, for which the old PWD works, is j = 2. Under this limitation 
we can compare the old and new PWD. We choose A = 5.8 m n in eq. ( fZBD and take as NN 
force the CD Bonn (np) |18| force restricted to j < 2, too. 
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In Fig. |3] we show two polarization observables, which show a big effect of the 3NF. (The 
differential cross section is only marginally affected by the 3NF.) Fig. [| demonstrates that 
the predictions for the old and new PWD do completely overlap. This is true for all other 
elastic observables, too. 

For the breakup process we calculated cross sections and all analyzing powers, spin 
correlation coefficients and vector spin transfer coefficients for standard kinematics. Again 
we chose as examples in Fig. |] observables with big 3NF effects. They are both observables 
under np QFS conditions. Again the agreement of the two curves for the old and new PWD 
is very good. This is generally true. If one chooses j max = 4 for the inner basis states 2 
and 3 instead of j max = 5 the deviations in the observables to the full result are up to a few 
percent. 

In order to show that also the extension of the new PWD to higher partial waves works 
well we performed a calculation of the triton binding energy with j max = 4. As NN force we 
used here the AV14 potential Our result for the triton binding energy is —9.36 MeV, 



which agrees exactly with the result of the Los Alamos - Iowa group [II 



We obtained the results presented up to now not by using eq. fl3"T| ) but the form 

.(pqalV^ (1 + P) \p'q'a\ 
= £ .(pqa] \p"q"a") x x {p"q"a"\ (1 + P) {p'q'a'), (94) 

This is identical to eq. (|3lD if the a" summation is unlimited. In the calculations, however, we 
restricted the double primed channels to the physical channels, which is an approximation. 
The old code for the 3NF did anyhow not allow the higher partial waves. Now using eq. (|31"D 
and the new PWD we can check the quality of that approximation. We compare the results 
for triton binding energies based on the new PWD of the 3NF achieved via eq. fl3i"D (denoted 
by E^ ] ) with the fully converged a"-sum to the result gained through eq. ( p4]) (denoted by 
Efr '), where the a"-sum is restricted to the same number of partial waves as used in the outer 
states. The latter number is fixed by choosing a certain maximal total two-body angular 
momentum j max - Our results are shown in table |. For j max = 5 we achieved convergence, 
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since the two binding energies E£ } and E ( b z > are equal. The effect of the truncation is large 
for j max = 2, whereas under the same restriction of outer partial waves is already close 

(2) 

to the final result. We also see, that the j max = 4 result for E b is close to the correct 
number E^\ 

The results with the NN force only are nearly converged at j max = 3, whereas the 
interplay with the 3NF causes a decrease in binding energy for j max = 3 and only for 
jmax = 4 the final result is essentially reached. 



V. SUMMARY AND OUTLOOK 

We demonstrated in section || that the way used up to now to decompose 3NF's in 
momentum space into partial waves leads unavoidably to numerical instabilities. They 
occur necessarily if the momenta and the angular momentum quantum numbers are high. 
The applications carried through up to now are not affected by that numerical instability, 
since the 3NF has been applied up to now only in low partial waves. For higher energies, 
however, high angular momenta are activated and the old PWD cannot be used. 

In section [TTT] we presented a new PWD. The basic idea is to treat the 3NF in the example 
of the 27r-exchange as if it would be a sequence of two pion exchanges among different 
pairs. This requires the choice of different Jacobi variables adopted to the two different two- 
body subsystems. Consequently in between there has to occur a recoupling of the Jacobi 
variables. In this manner we not only achieve a much simpler analytical form of the partial 
wave decomposed 3NF, but also an expression, which is free of any artificial singularity 
and numerically stable in all angular momentum states. The recoupling between the two 
meson exchanges and also to the external basis states, which require the third possibility 
for the choice of a two-body subsystem, are performed with the aid of certain permutation 
operators, which are laid out in appendix A. The application of these permutation operators 
lead to two-dimensional interpolations, which require very efficient and accurate tools. To 
that purpose we have rewritten standard hermitian cubic splines into a new form, which is 
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shown in appendix B. 

Finally we have shown in section [IV] two numerical examples for the reliability of the 
new PWD of the 3NF. We used the TM 3NF model and calculated in the first example 
observables in elastic nd scattering and the breakup process at Ei a b = 10.3 MeV. The forces 
were restricted to j max = 2 so that the old PWD could also be used. For all observables we 
found a very good agreement between the predictions achieved with the old and new PWD. 
In the second example we extended the new PWD of the TM 3NF to higher partial waves 
with j max = 5 and calculated the triton binding energy. We could confirm the result of the 
Los Alamos - Iowa group truncated at j max = 4. The convergence of the triton binding 
energies as a function of j max as displayed in table [I] shows, that in contrast to the case of 
NN forces only, the interplay with the 3NF can lead to a more oscillating approach of the 
limiting value. More benchmark calculations are underway. 

Summarizing we can say that we are now able to calculate the 3NF up to any partial 
wave, limited only by computer resources. This new scheme can of course be taken over 
from the 27r-exchange, exemplified here, to any other two-meson exchange. 

It is interesting to note that the new PWD has still another application. If one wants 
to calculate three-body meson exchange currents as they occur for example in electron 
scattering on 3 He ( 3 H) or heavier nuclei one faces exactly the same problems as described 
here. In order to get correct results for all partial waves one has to follow for the PWD of 
the three-body meson exchange currents the scheme which we have introduced above. 
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APPENDIX A: THE PERMUTATION OPERATORS 



In the Faddeev scheme [12] [I3| we apply the permutation operator P = P12-P23 + -P13-P23 
always to states which are anti-symmetric in the two-body subcluster 2-3. Regarding 
eqs. ( |55D and (0) we see the recoupling matrix elements x {pqa \ p'q'a') 2 , 2 (p<l a J I v'l' a 'j)^ 
and 3 (pqa \ p'q'a') v Since |) 3 = P13-P23 |)i an d |) 2 = P12-P23 |)i one easily finds that all three 
matrix elements are of the type t (\ P12-P23 The expressions given in fl2[] and |§ for P 
are in fact evaluated there with the help of 1 (| P12-P23 |)i and an originally occurring phase 
factor (— y'+ s '+ t ' has been replaced by (—1), since there only antisymmetric two-body states 
have been used. Redoing that replacement one finds the desired matrix element from the 
one of P given in |12j and 0| by multiplying it with |(— y'+ s '+ t '+ 1 . For the states without 
iso-spin one has just to drop the iso-spin factor. 



Out of the many various forms to evaluate permutation operators |2(| we need here three 
of them. The recoupling between the quasi two-body force quantities requires 



x {pqaj\ P 12 P 23 \p q'a' J ) l = - (-)< +s +1 j ^ dx V ^ +2 p/l , +2 ( 



J a ., 



with 



7T = W -q 2 + q' 2 + qq'x 



(Al) 



7r' = Jq 2 + -q' 2 + qq'x 



(A2) 



and 



G aja >(qq'x)=j:Pk(x) E E 1 h+l ' 2 q' h ^9 k hhl ' A 



h+h=l l[+l' 2 =l' 



The purely geometrical quantity g^ 1 ^ 2 / 1 ' 2 has to be taken now without iso-spin: 
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(A3) 
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x k 
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\ (2/0! (2J 2 )I \ (2JJ)! (2Z 2 ) 



(2Z' + 1)! 



r /i z 2 / 1 f z 2 z; /' 

//' i A L f ) I A' L f 
f k L 

x < } C(fcO ZiO, f'0) C(k0 Z 2 0, /0) 

/' Z 2 fc 



(A4) 



The very left recoupling coefficient has to be taken in the form where the 5-functions act 
both to the right: 



i 



1 (pqa\P 12 P 23 \p'q'a') 1 = - J ^ ,/,- 



G aa '(pqx) (A5) 



with 



* = V 4 P + f6 g + 4^ 



x = \jp 2 + -^q 2 -pqx 



(A6) 



and G aa /(pqx) can be taken for instance from ||. 

Finally the recoupling coefficient from the states of type 1 to 3 require a form where the 
^-functions act both to the left: 

(pqa\ P 12 P 23 \ P 'qW) = \ dx %=^> G^p'q'x) (A7) 



<pl-\-1 q\-\-1 



with 



1 „ 9 , 3 
^ 16 g 



X = ]lp' 2 + -q' 2 + p'g'x 



(A8) 



In eq. ( |A7f) we need the same quantity G as in eq. (|A5[) , but with interchanged arguments. 



APPENDIX B: CUBIC HERMITIAN SPLINES 



In appendix A we encountered two-fold interpolations. They have to be performed 
for very many channels and should be as optimal as possible. We found that our usual 
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basis splines JTjj are not efficient enough to perform these two-fold interpolations within a 
reasonable time. Therefore we need a more efficient interpolation algorithm. 

Such an algorithm can be found by using cubic hermitian splines. Although cubic her- 
mitian splines are well known in the literature (see for example ||21||), we will give here a 
short introduction in order to explain our way to use them. 

Consider a function f(x) given at certain grid points x», i — 1, . . . , n. Let x be positioned 
in the interval [xj,Xj+i]. For the sake of simpler notation we call the end points Xj = x\ and 
Xj + i = x 2 . Then one defines an unique cubic polynomial fi(x) by the following constraints: 

fi(xi) = /Oi) 
fi(x 2 ) = f(x 2 ) 
fl( Xl ) = f(x 1 ) 

fl(x 2 ) = f'(x 2 ) (Bl) 

Therefore these interpolating functions fi(x) and their derivatives f[(x) are continuous at 
the grid points Xj. They are given by 

fi(x) = fix^Mx) + f(x 2 )Mx) + f'MMz) + f'MMz) (B2) 
in terms of the spline functions 

<pi{x) = ( ,^ 2 _ ^ ((x 2 - Xi) + 2(x - Zi)) J 

\ ^2 X\ J ) 
MX) = ^ ~ %l) + ^ ~ 

(x- x 1 )(x 2 -xf 

h{x) = 



(x 2 - Xi) 2 
(x — x\) 2 (x — x 2 ) 
(x 2 - Xif 



Mx) = — TTi ( B3 ) 



We approximate the derivatives f'(%i) and f'(x 2 ) with the help of a quadratic polynomial, 
which is uniquely defined by the function values at the grid point and its two neighbours. 
Calling x i _ 1 = x and x i+1 = x 3 we get 
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/ (Xi) « /(x 2 ) /(x ) 



+ /(Xl) 



X 2 — X\ X\ — Xq\ 1 



X\ — Xq X2 — X\ / X2 — Xq 



f{x 2 ) ~ f(x 3 ) - J /(xx) -J— 

X3 — X2 3? 3 — X\ X2 — X\ X3 — X\ 

+ /( x 2 ) ( 512 - *z£i) _J_ ( B4) 

\X 2 -Xi X3-X2/ X3 — X1 

At the end points Xi and x n we define the quadratic polynomial by /(xi), /(x 2 ) and /(x 3 ) 
and by /(x n _ 2 ), /(x n _x) and /(x n ), respectively This is achieved by putting x = x 3 in the 
first case and X3 = x n _ 2 in the second case. 
Insertion of eq. (TBI) into eq. flB2|) yields 



A(i) = E^) /(•'•,) (B5) 
i=o 

for the interpolating function in the z'-th interval of the grid points. Thereby we are lead to 
the modified spline functions 

x 2 — X\ 1 



S Q (x) = -<p 3 (x) 



X\ — Xo X 2 — Xq 

'x 2 — X\ Xi — X \ 1 / v x 3 — x 2 



S\(x) = <j)%{x) + fa(x) (— — — -— — —\ 4 (x) 

VXi — Xn X 2 — Xi / X 2 — X 



v Xi — Xo X 2 — X\J X 2 — Xo X 2 — Xi X3 — Xi 

a t \ x 1 \ 1 x 1 \ Xl ~ x ° ^ x r \ f x 3 ~ X2 x 2 — Xi\ 1 
S 2 (x) = fa{x) + </> 3 (x) h04W 

X 2 — Xi X 2 — Xo \X 2 — Xi X3 — X 2 / X3 — Xi 

S 3 (x) = 4 (x) — L_ (B6) 

X3 — X 2 X3 — Xi 

Eq. ( |B5|) is very well suited for the numerical usage. The modified spline functions Sj(x) 
are independent of the function values and depend only on the grid points and the actual 
position x. Therefore they can be prepared beforehand. This is very important if one has 
to interpolate very many functions given at the same grid points as we have to do in our 
3NF code. 

The form of eq. ( |B5| ) is the same as the one found in |14j for basis splines. The difference 



is, that the sum for the basis splines runs over the whole grid, whereas the sum for the 
hermitian splines in eq. ( [B5| ) has only four terms related to the four grid points nearest 
to the interpolation point x. (Basis splines are global splines, whereas hermitian splines 
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are local.) Assuming a grid of typically 20 points the one-dimensional interpolation using 
hermitian splines needs only | operations compared to basis splines. For a two-dimensional 
interpolation the number of operations is reduced by a factor of ^. 

For the two-dimensional interpolation one has to make a bi-cubic ansatz for the inter- 
polating functions fij(x,y). To define a bi-cubic function uniquely we need 16 constrains, 
which we choose as 

fij(xi,yi) = 



9fij(x,yi) 



dx 
dfij(x 1: y) 



df(x,yi) 



dy 

d 2 fij{x,y) 



X=X\ 



y=yi 



dx 
df(xi,y) 



dx dy 



dy 
d 2 f(x,y) 



X=Xl 



y=yi 



x=x lt y=yi 



dx dy 



(B7) 



x=x 1 ,y=y 1 



and identical expressions for the other three points (xi, 1/2), (x 2 , yi), and (x 2 , 1/2), respectively. 
Hereby the four grid points (xi,yi), (xi,y2), (x2,yi) and (#2,2/2) are nearest neighbours to 
the interpolation point (x, y) in the xy-plane. 

The partial derivatives are approximated as in the one-dimensional case. The second 
derivative is estimated by a bi-quadratic polynomial, which is uniquely given by the function 
value at the specific grid point and the function values of the eight surrounding grid points. 

Following these steps one yields 



k=0 1=0 



(B8) 



for the interpolating function. It is a fairly easy exercise to show that the two-dimensional 
spline functions Sff(x,y) are simply given by 



s£\x,y) = S M {x)S l (y) 



(B9) 



Analogical equations hold for interpolations in more than two dimensions. 

To our experience one- and two-dimensional interpolations using hermitian splines are 
at least of the same accuracy as interpolations based on basis splines. 
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TABLES 



3max 


tpNN 

h b 


4" 


E f> 


2 


7.58 


9.34 


9.51 


3 


7.67 


9.27 


9.25 


4 


7.68 


9.34 


9.36 


5 


7.69 


9.32(9) 


9.32(4) 



TABLE I. Triton binding energies calculated using the AV14 NN force alone (Eg N ) and to- 
gether with the TM 3NF (A = 5.8m w ). For explanation of E^ and Ef ] see text. 
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FIGURES 
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FIG. 1. The 2vr-exchange 3NF. 
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FIG. 2. The part F 4 (1) of the 27r-exchange 3NF. 
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FIG. 3. The spin correlation coefficient C yy and the vector to tensor spin transfer coefficient 

Ky' z ' for elastic scattering at Ei a b = 10.3 MeV. The solid line is the result without 3NF, the 

overlapping long and short dashed lines are the results with 3NF using the old and new PWD, 

respectively. 
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FIG. 4. The differential cross section and the deuteron analyzing power A y under np QFS 

conditions for the breakup process at Eiab = 10.3 MeV as functions of the scattering angle of one 

nucleon. Description as in Fig. 3. 
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